Variational Monte Carlo Study on the Superconductivity in 
the Two-Dimensional Hubbard Model 

Kunihiko Yamaji'^'^I, Takashi Yanagisawa", Takeshi Nakanishi^H and Soh 

Koike*'" 

" Electrotechnical Laboratory, 1-1-4 Umezono, Tsukuha, Ibaraki 305-8568, Japan 
^ Institute of Materials Science, University of Tsukuba, Ten-nodai, Tsukuba 305-8573, Japan 

(Received May 29, 1998) 

The possibihty of superconductivity (SC) in the ground state of the two-dimensional 
(2D) Hubbard model was investigated by means of the variational Monte Carlo method. 
The energy gain of the d-wave SC state, obtained as the difference of the minimum 
energy with a finite gap and that with zero gap, was examined with respect to depen- 
dences on U, electron density p and next nearest neighbor transfer t' mainly on the 
10 X 10 lattice. It was found to be maximized around U = 8 (the energy unit is nearest 
neighbor transfer t). It was shown to sharply increase for negative values of t' and 
have a broad peak for t' ~ —0.10. For these value of t' the energy gain was a smooth 
increasing function of p almost independent of the shell structure in the region starting 
from ~ 0.76 up to the upper bound of investigation 0.92. This clearly indicates that 
the result is already close to the value in the bulk limit. For t' = 0, the energy gain de- 
pended on the electronic shell state. This suggests the 10 x 10 lattice is not sufficiently 
large for this case, although it is highly plausible that the bulk limit value is finite. 
Competition between the SC and the commensurate SDW states was also investigated. 
When t' = 0, the ground state is SDW in the range of p >~ 0.84. The SC region 
slightly extends up to ~ 0.87 for t' ^ —0.10. Consequently the present results strongly 
support an assertion that the 2D Hubbard model with t' ^ —0.1 drives SC by itself in 
the p region from ~ 0.76 to ^ 0.87. The above features are in a fair agreement with 
the phase diagram of the optimally and overly hole-doped cuprates. The energy gain 
in the SC state with suitable parameters is found to be in reasonable agreement with 
the condensation energy in the SC state of YBa2Cu307. The corresponding t-J model 
proves to give an order-of-magnitude larger energy gain, which questions its validity. 
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§1. Introduction 

Recently the mechanisms of superconductivity (SC) in high-temperature cuprate superconductors 
and organic superconductors have been extensively studied using various two-dimensional (2D) 
models of electronic interactions. The 2D Hubbard model is the simplest and the most fundamental 
one among such models. Early studies of this model using the quantum Monte Carlo (M.C.) method 
indicated the existence of an attractive interaction for SC |T[]. However, this possibility has been 
controversial. Some authors asserted from quantum M.C. results that the enhanced SC correlation 
does not develop into the predominant one at low temperatures or in the ground state of this model 
[|,||,§,|5|. Some authors supported this possibility by numerical results using the variational Monte 
Carlo (M.C.) method |^. Since it is of prime importance to establish the simplest electronic model 
for superconductivity and since these early studies had shortcomings due to the restrictions to 
the investigated parameter space, the problem deserves more extensive investigation from various 
aspects such as dependences on the system size, the on-site Coulomb energy U, the electron density 
p and the next nearest neighbor transfer energy t'. Recently, appropriate values of t' were found 
to remarkably enhance SC correlations 0, ^0]. We apply the variational Monte Carlo method 
mainly to the 10 x 10-site system with electron numbers from 68 to 92. This method has 
a merit that it allows to treat larger values of U than the quantum M.C. methods. We revisit the 
issue investigated in ref. H, carrying out more extensive study with higher precision in a wider 



parameter space. In a previous preliminary report |T2[ we worked out a variational Monte Carlo 
calculation 84 electrons on the 10 x 10 lattice and confirmed that the 2D Hubbard model has the 
d-wave SC state. In the present work first we examine the [/-dependence of the energy gain due to 
the condensation into the SC state and find the optimal value of U is about 8 in units of the nearest 
neighbor transfer energy t. In the previous report we confirmed that the energy gain in the SC 
state sharply increases with the introduction of t' with a negative value. In this report we calculate 
the energy gain as a function of t' in a sufficiently wide range for a fixed value of f/ = 8 and show 
that it has a broad peak around a negative value t' ~ —0.1. The approximate peak energy gain in 
the SC state, i.e., at t' ~ —0.1, is plotted as a function of electron density p together with that for 
t' = 0. When t' ~ —0.1, the energy gain starts to be finite at about p = 0.68 and increases fairly 
smoothly with increase of p in the range of 0.76 < p < 0.92. We carried out calculations in both 
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cases of open and closed shells. The smoothness of the energy gain against p indicates that the 
dependence of the energy gain on the shell structure is weak, in contrast to the case with t' = 0, 
and that the results are already close to the bulk limit. This means that the 2D Hubbard model 
with t' ~ —0.1 has a definite bulk-limit energy gain per site in the SC state. While when t' = 0, 
the energy gain in the closed shell state was found much smaller than that in the open shell state. 
This indicates that the 10 x 10 lattice is not sufficiently large for the case of t' = 0, contrary to the 
case of t' ~ —0.1, although a kind of average of the results in both shell states is very probable to 
give a finite energy gain per site in the bulk limit. Of course it is needed to treat larger lattices 
with t' = but some important findings at the present stage are considered to deserve a report, 
with the latter problem left to be settled, since with larger sizes computation needs extremely long 
time. 

Thirdly, we calculate the energy gain in another ordered state, i.e., spin density wave (SDW) 
state. We treat the SDW with a fixed wave vector {tt,tt). Here the length unit is the lattice 
constant. The energy gain in the SDW state was found to quickly drop with decrease of p from 
unity and to vanish at p just below 0.84 when t' = 0. However, at p = 0.84 the SDW state was 
slightly more stable than the SC state so that the boundary between both states lies just below 
p = 0.84 for t' = 0. t' is known to destabilize the SDW in the hole-doped state, while it promotes 
SC pairing. We study to what extent the phase boundary between the competing SC and SDW 
states is affected by the negative t'. The boundary is appreciably pushed up to the higher p side 
with increase of the absolute value of the negative t', although not largely. 

We point out that the resultant SC region from p ~ 0.76 to p ~ 0.86 and the increased SC energy 
gain with increase of p are in accord with experimental features of cuprate high-Tc superconductors 
in the overdoped region. The observed energy gain in the SC state is found to be in fair agreement 
with the experimental condensation energy estimated from the critical magnetic field He and the 
specific heat of YBa2Cu307. On the other hand the t-J model as an effective Hamiltonian of the 
Hubbard model is noticed to give an enormous overestimation of the condensation energy. 

In the next section the model and the method are described. Results on the SC and SDW ground 
states are given in §3 and §4, respectively. In §5 the obtained results are compared with experiments 
and also with other computational works. The final section gives the summary. Short reports on 
the present results were given in [|12|,[1^. 

§2. Model and Method 

Our model is the 2D Hubbard model defined by 

H = -t + H.c.) +UJ2 cltC^TcJiCii' (1) 

<jl>,o- j 



3 



where (cjo-) is the creation (annihilation) operator of an electron with spin a at the jth site; the 
sites form a rectangular lattice; t is the transfer energy between the nearest-neighbor (n.n.) sites; 
t is our energy unit; < j7 > denotes summation over all the n.n. bonds. V is the on-site Coulomb 
energy. In this report, we also study the effect of t' between n.n.n. sites by including 

Hnnn = -t' ^ (cJ-Qa + H.C.) (2) 
«jl»,cr 

in the Hamiltonian; in the above equation « jl » means summation over the n.n.n. pairs. 
Our trial wave function is a Gutzwiller-projected BCS-type wave function defined as [p!0|,p 



= Pg^BCS, (3) 
k 

where Pq is the Gutzwiller projection operator given by 

^G = nii-(i-^?KT^d; (5) 

i 

(yf is a variational parameter in the range from to unity and j labels a site in the real space. 
P/Vg is a projection operator which extracts only the states with a fixed total electron number Nf.. 
Coefficients and appear in our calculation only in the ratio defined by 



= — 2t(cos kx + COS ky) — At' cos k^ cos ky — /i, (7) 

where is a fc-dependent gap function defined later; /i is a variational parameter working like 
the chemical potential in the trial wave function; c^^ is the Fourier transform of Cjv. Neglecting 
constant factors, \E's can be rewritten as 

VI/, ~ Piv.PGexpE(^fc/«fe)4t^^fc|]|0 >' 
k 

= P^,PGexpEa(j,/)44]|0 >, (9) 
- PG[T.^ijJ)c]Af'^"\^>^ (10) 

= PG Y1 Mjl^--JNe/2-Jl,--jNe/2) 

Jli--,ijVe/2:'lv,'jVe/2 

X cj tctf-cj tC/" |C;t I •••cj' , ||0>, (11) 



where a(j, /) is defined by 



a(j, /) = (1/iV,) exp[ife ■ {Ri - Rj)], (12) 

k 
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with Ng being the number of sites and 



^(jl, J2, ■■■,jN,/2', h, h, = 

«(j2,/l) a{j2,l2) ■■■ a{j2,lNj2) 



aUN,/2,h) a{jN,/2j2 



Ne/2, lNe/2 



(13) 



ji, 72, • • • are arranged in the ascending order; so are /i, /2, . . .. Then the ground state energy 



Eg =< H >=< ^s\H\^s > I < ^sl^s > 



(14) 



is obtained using a M.C procedure |TO|,|Il|. In order to minimize computation time in the M.C 
computation the values of the cofactors of the matrix in eq. (13) were stored and corrected at 
each time when the electron distribution was modified. Errors accumulated in the cofactors after 
corrected many times were avoided by recalculating cofactors as determinants, instead of getting 
them by correction, after a certain number of M.C. steps. We tested our programs using exact 
diagonalization results for small systems. We optimized Eg with respect to g, and fi. In the 
later stage of work we employed the correlated measurements method |1^,|T5| in the process of 
searching optimal parameter values minimizing Eg. It enables to precisely locate the minimum 
position. 

§3. Results for the Superconducting Ground State 

3.1 Case of Simple 2D Hubbard Model 
We studied the cases of the d-, extended s- (s*-) and s-wave gap functions as follows: 



d Af^ = A(cos A;^. — cos fc. 



A (cos kx + cos ky), 
A 



(15) 



s Af^ = 

The sizes of the lattice we treated are 6x6 and 10 x 10 having electron density close to unity with 
slight hole doping into the half-filled state. 

Results indicating the occurrence of the (i-wave superconductivity were obtained even for the case 
of Ne = 32 on the 6x6 lattice with the periodic and the antiperiodic boundary conditions (b.c.'s) 
for the X- and the y-direction, respectively, with U = 8 and t' = 0. This set of b.c.'s was chosen so 
that A^ does not vanish for any fe-points possibly occupied by electrons; otherwise zero division 
occurs in eq. (0). Eg/Ng is minimized at A ~ 0.10. Here g = 0.3038 and fi = —0.48. The energy 
gain due to the SC gap formation, i.e., SC condensation energy, was estimated at ~ 0.00028/site 
from the difference between the minimum and the intercept of the Eg-veisns-A curve with the 
vertical axis. 
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Fig. 1. Computed ground state energy per site Eg/Ng is plotted against A for the case of 84 electrons on 
the 10 X 10 lattice with C/ = 8 and t' = 0. Filled circles are for the d-wave gap function with g and fi 
optimized for each A. Filled squares and triangles are for the s*- and s-wave gap functions, respectively. 
The diamond shows the normal state value. Straight lines between data points are the guide for the eye. 
The thin curve is a parabola given by the mean-squares fit to the d-wave data. 



Main results written in this report were obtained for the 10 x 10 lattice with the same boundary 
conditions. Calculated energies per site for several types of wave functions with A^e = 84 on the 
10 X 10 lattice are shown in Fig. 1 for the case of C/ = 8 and t' — Here Eg/Ng is plotted as a 
function of A for the three types of gap functions given in eq. (15). With the lattice, the b.c.'s 
and A^e = 84, the electronic shell structure in the hmit of C/ = is open, i.e., some fe-points are 
partially filled, as is illustrated in Fig. 2, which displays the electron occupancy at the fc-points. At 
each value of A in Fig. 1, g — 0.30 was chosen as the initial value of g and then the optimal value 
of // was found by the least squares fit of Eg as a function of // to a parabola. Using this value of 
/X, g was optimized again. Since Eg was a smooth function of g, the obtained optimal g and jJL are 
sufficiently accurate. Using these parameter values. Eg was obtained as the average of the results 
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Fig. 2. Distribution of fe-points in the reciprocal lattice in the case of the 10 x 10 lattice with the boundary 
conditions mentioned in the text. Closed and empty circles denote doubly occupied and empty sites, 
respectively, and circles partially closed are partially occupied sites in the case of 84 electrons. A, B, G 
label fe-points lying in the neighborhood of the Fermi energy for later use. 



of eight M.C. calculations each with 5 x 10^ steps at A = 0.01,0.04 and 0.08 for the rf-wave. The 
standard deviations were less than 0.00011. At other points, the numbers of M.C. calculations and 
steps were different but their error bars were within 0.00015 with the total M.C. steps greater than 
2 X 10^. The diamond shows the normal state value, —0.73585 ± 0.00024, obtained from 20 M.C. 
calculations each with 10^ steps, using the Gutzwiller-projected normal state wave function. In the 
employed normal state (vr/G, 77r/10) and (— 7r/6, — 77r/10) are fully occupied but (vr/G, — 77r/10) and 
(— 7r/6, 77r/10) are unoccupied in the component of the trial wave function prior to the Gutzwiller 
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projection. 

Clearly, Eg/Ng is minimum at a finite value of A ^ 0.08 in the case of the ci-wave gap parameter. 
The optimal parameter values at A = 0.08 are g = 0.3037 and /i = —0.4263. The least squares fit of 
a parabola to Eg/Ng as a function of A is of good quality, as seen in Fig. 1 and gives the minimum 
at A = 0.082. The curves of Eg/Ng for the s- and s*-wave gap functions have definite positive 
slopes at small A and are extrapolated for A = to ~ —0.7354 and ~ —0.7353, respectively, 
which are practically equal. These values are slightly higher than the value ~ —0.73605 given by 
extrapolating the ci-wave fitting parabola to A = 0. The normal state value of Eg/Ng = —0.73585 
lies between the two groups of extrapolated values. The differences are due to a kind of size effect, 
as was explained in the previous report [|l^. These differences in the case of 10 x 10 lattice are 
much smaller than the depth of the minimum of the d-wave curve. 

The energy gain per site in the rf-wave state was obtained from the parabola fitting to the energy- 
versus-A curve as its depth of the minimum in reference with the intersection of the fitting curve 
with the vertical axis as ~ 0.00155/site. It is five times larger than that for A^e = 32 and A^^e = 36. 
This suggests that the size effect is considerable in the 6x6 system. This energy gain per site 
in the case of Ng = 100 and A^^e = 84 is of the order of magnitude of the BCS superconducting 
condensation energy ~ (state density) xA2 ^ 0.0010 since the state density per site per spin around 
this density is approximately equal to l/27r in units of t. 

In order to check the superconducting nature of the ground state with a finite value of A, the 
correlation functions of BCS pair operators were calculated. Superconducting pair correlation 
functions Dap^l), a, P = x, y, are defined as: 

D^p{l) =<Al{t + l,j)Ap{t,j) >, (16) 

where Aa{i,j),(y = x,y, denote the annihilation operators of singlet electron pairs staying on n.n. 
sites as: 

^x{i,j) = QjiQ+iji - dj^d+iji, (17) 
^yihj) = CiiiQj+n - Cij|Qj+i|, (18) 

where Cjjo- means the annihilation operator at site with spin cr. The average < . . . > is defined 



in eq. (0). The result for a state close to the minimum Eg point in Fig. 1, i.e., A = 0.078, 
fi = —0.428 and g = 0.30, is shown in Fig. 3. The vertical scale is enlarged to highlight the long 
range parts. The correlation extends over the whole lattice as expected, showing a clear contrast 
to the normal state. The ci-wave nature appears in the negative sign of Dxy{l) for / = 2 ~ 5. 

3.2 U -Dependence 

The [/-dependence of the energy gain in the SC state was shown in Fig. 4. with the same values 
of Ng and N^. The energy gain was determined from the depth of the minimum, as explained 
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Fig. 3. Parts of the correlation functions Da/3 of superconducting pair operators for case of the d-wave 
minimum state shown in Fig. 1. The horizontal axis is the distance I between two positions of two pair 
operators. Closed symbols are for the superconducting state and open ones for the normal state. Dap is 
defined by eq. ([l6|). 



above. Except for the above-mentioned case of f/ = 8, the minimum position defined by optimal 
A, g and /i was determined by means of the correlated measurements and Eg was calculated for this 
position with the same number of M. C. steps as for U = 8. Next, Eg for A = 0.01 was obtained 
and then the fitting parabola allowed to get the energy gain. The energy gain in this definition 
gives the condensation energy in the SC state in the present approximation. The energy gain has 
the maximum for about U = 8. It quickly decreases with decrease of U. It gradually decrease with 
increase of U over 12. In contrast, the value of optimal A is nearly proportional to t/ up to f/ = 
12 and then saturates for larger U. 

3.3 t' -Dependence 

In the 2D Hubbard and the 2D d-p models, the one-electron state density increases peak-wise 
around the energy of the van Hove singularities located in the fe-space region around (0, vr) and 
(vr, 0). When t' is negative in the 2D Hubbard model with Hnnn the energy level of van Hove 
singularity moves toward the Fermi energy in the hole-doped systems, as is seen from eq. (7), which 
realizes a situation favorable to get superconductivity for most theories. In fact Shimahara and 
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Fig. 4. Energy gain Ai?g (closed circles) and gap parameter A (open circles) in the superconducting state 
are plotted against on-site Coulomb energy U for the same system as in Fig. 1. 



Takada showed that in the RPA framework Tc increases with the introduction of t' |T^ . Importance 



of t' in the physics of high-Tc cuprates was pointed out in • From the viewpoint of the two-band 
mechanism of superconductivity in which the pair-wise transfer of electrons from a certain fc-region 
to another fe-region promotes SC, the location of the fc-space region highly contributing to the 
state density in the neighborhood of the Fermi energy should lead to a larger energy gain in the SC 
state [|18| , p!9| . The state density peak is known to be further enhanced with the increase of electron 
correlation, as the single-particle dispersion along the lines from these points to (0, 0) becomes 



anomalously weak pO|. The increased state density around the singularity has been argued to 
enhance Tc of the (i-wave superconductivity [^. The value of t' was argued to correlate with the 
maximum Tc values of cuprate subfamilies Quantum M. C. studies indicated that i! enables 

the bulk superconductivity in the 2D Hubbard model 0. 

We have examined the t'-dependence of the energy gain, taking account of Hnnn in (0) in our 
model. In the case of A^g = 32 on the 6x6 lattice, the energy minimum became slightly deeper 
with the change of t' from zero to —0.25, but at t' = 0.25 the minimum became shallower by a 
factor of 2. In the case of Ne = 84 on the 10 x 10 lattice the energy gain was maximized around 
t' = —0.1 ~ —0.15 as shown in Fig. 5. The error bar for AEg/Ns in the figure is about 0.0003. 
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Fig. 5. Energy gain per site AEg/Ng in the superconducting state is plotted by closed circles with the left 
vertical scale as a function of next nearest neighbor transfer energy t' in the case of 84 electrons on the 
10 X 10 lattice with {7 = 8. Values of A are shown by open triangles with the right-hand-side vertical scale. 



A is also displayed in the figure with the right-hand-side vertical scale. The error bar for A is 
about 10 percents of the magnitude. With a positive value of t' — 0.10, the energy gain quickly 
decreases, nearly vanishing when t' — 0.25. It also nearly vanishes for t' ~ —0.40. Figure 6 shows 
the SC pair correlation functions as functions of the distance between pair operators in the optimal 
state for the case of t' — —0.10. At the farthest distance their values are twice larger that those for 
t' — (Fig 3), as expected. Thus, an appropriate negative value of t' enhances the d-wave SC in the 
hole-doped case in the Hubbard model, which is in qualitative agreement with the above-mentioned 
references. The i'-dependences of the energy gain and A in the case of — 80 exhibited in Fig. 
7 show similar features as in Fig. 5. 

When one plots one-electron levels (Fig. 8), £^ = + /x, for fe-points in the non-interacting 
limit as a function oit', a few occupied levels and a few empty levels are found to be bunched into a 
close neighborhood of the Fermi energy sp within ±0.07 as t' takes —0.1 ~ —0.15. These levels are 
located in the fc-space mainly in the neighborhood of the level of the van Hove singularity. Since each 
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Fig. 6. Parts of the correlation functions of superconducting pair operators for the case of the largest SC 
energy gain in Fig. 5 with t' = -0.10, A = 0.0941, g = 0.3003, and fi = -0.6152. The notation is the 
same as in Fig. 3. 



one-electron level is orbitally four- (or rarely two-) fold degenerate due to symmetry, the number of 
one-electron states whose levels are located in the neighborhood of ep is considerable. Since such 
a feature is known to enhance the SC pair correlation functions in the two-chain Hubbard model 



1^ even in strongly correlated situations, this feature is considered to bring about the remarkable 
increase of the energy gain. Although the particular peaking location of t' may be slightly size- 
dependent since the one-electron levels are size-dependent, such a concentration of one-electron 
levels around ep is considered to occur for any system size when the appropriate negative t' pushes 
down the van Hove singularity close to ep. Therefore, this effect should work even in the bulk limit. 

3.4 Electron- Density Dependence 

Figure 9 shows the energy gain per site in the SC state as a function of the electron density 
p = Ne/Ng for the cases of t' = 0, —0.10 and —0.15. The corresponding values of A as a function of 
p are exhibited in Fig. 10. When t' = —0.10 whose curve is exhibited with closed circles, the cases 
of p = 0.68, 0.76, 0.84 and 0.92 have an open shell, in another word, there are partially filled k- 
points in the electronic state in the U = limit, as is shown in Fig. 2 for A^s = 10 x 10 and = 84 
with t' = 0. While the other cases with p = 0.72, 0.80 and 0.88 have a closed shell. The curve 
has almost no dependence on the shell state in the region of 0.76 < p < 0.92, making a definite 
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Fig. 7. AEg and A as functions of t' in the case of p = 0.80. Other parameter values are the same as in 
Fig. 5. 

contrast with the case with f = described later. This suggests that the curve is already close to 
the bulk limit one in this region. For t' = —0.15 with open diamonds, the cases of p = 0.68, 0.76 
and 0.88 have an open shell while others have a closed shell. Here we observe a weak dependence 
on the shell structure. However the deviation from the rough average between the open-shell and 
closed-shell curves are only of the order of 1/10 of the average energy gain per site in the region of 
0.76 < p < 0.92, which is small enough to judge that the points are already close to the bulk limit. 
The smooth dependence on p is also observed in the SC gap amplitude A for t' = —0.10 and —0.15 
in Fig. 10. These results clearly support a statement that the 2D Hubbard model with t' ~ —0.1 
gives a finite condensation energy for the rf-wave SC state in the bulk limit in the region of p from 
~ 0.76 to ~ 0.86. For p = 0.72 a minimum with finite A was obtained in both cases of t' ~ —0.1. 
But the SC energy gain was within the error bound so that AEg/Ns ~ 0. For p = 0.68 a finite 
value was obtained for both cases but it is so small that it may not remain finite in the bulk limit. 

Open circles for t' = in Fig. 9 form a zigzag curve as a function of p. Relatively larger are the 
energy gains for p = 0.76,0.84 and 0.92, for which the electron shell is open. The values for the 
cases of p = 0.72, 0.80 and 0.88 with closed shell are smaller. This shell dependence is a system-size 
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Fig. 8. One-electron levels of fe-points A, B, G in the neighborhood of ep are plotted as functions 
of t' . Due to orbital and spin degeneracies, each levels can accommodate eight electrons except the A level 
which can accept four. The dotted curve exhibits ep in the case of the U — Q limit. 



effect and indicates that the 10 x 10 lattice is not sufficiently large when t' = 0. It suggests also 
that the larger the SC energy gain is, the smaller the sufficient system size is. However, the average 
curve between the two curves formed of the points of both kinds of shell states is expected to be 
a highly probable approximation to the bulk limit and the model with t' = is also probable to 
drive SC, although more extensive calculation with large lattices is needed to remove the room for 
doubt. Incidentally, for p = 0.68, no finite energy gain was obtained although the shell is open. 

In Fig. 1 of the precursive work of Giamarchi et al. we find a corresponding value for the 
SC energy gain per site. It is for p = 0.8125 on the 8x8 lattice with U = 10 and t' = 0. Our 
value for p = 0.80 on the 10 x 10 lattice with U = 8 and t' = should be close to the value. 
The shell structure for p = 0.8125, or A^e = 52, on the 8x8 lattice is closed; so is the case with 
p = 0.80 on the 10 x 10 lattice is closed. However, our value is 3.38 x 10"'^, being about 1/3 of 
the above-mentioned value, and if we put it in their figure it is located in the lower end part of the 
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0.95 



Fig. 9. Energy gain per site AEg/Ng in the superconducting state against electron density p. Open circles 
are for t' = 0, closed ones for t' = —0.1 and closed diamonds for t' = —0.15. The lattice is 10 x 10 and 
U = 8. 



large error bar attached to the data point. In the above-mentioned figure the energy gain vanishes 
around p ~ 0.56 but in Fig. 9 of the present work it vanishes around p ~ 0.72. We beheve that 
the present results are more improved than the results of Giamarchi et al. because of improved 
accuracy of computation. 

§4. Competition with the SDW State 

4-1 Magnetic State of the 2D Hubbard Model 

When the density of doped holes is small or zero, the 2D Hubbard model takes an antiferro- 
magnetic state as its ground state. With increase of doped hole density, the magnetic order is 
destroyed and SC appears. At what hole density does the SC state start? We have investigated 
the transition between the pure SC and the pure uniform SDW states by computing the energy of 
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Fig. 10. Amplitude of the gap function A is plotted as functions of p for t' = 0, —0.10 and —0.15, corre- 
sponding to the AEg/Ns versus p curves in the preceding figure. The error bar is about 0.01 for p ^ 0.92. 
For p - 0.7 it is about 0.005. 



the SDW state by means of the variational Monte Carlo method. The trial SDW wave function is 
written as [|2|,[n| 

*SDW = ^G^SDW. (19) 
k k' 



Ui 



fc = [(1 - ^fc/V^fc + M^)/2Y'\ (21) 
^fc = [(l+«^fc/\A4+^)/2]'/', (22) 

^k = i^k-^k+q)/'^^ (23) 

^fc = ^fc + /^' (24) 

where Pq is the Gutzwiller projection operator defined by eq. (5). Summation over k and k' in 
eq. (20) is performed over the filled fc-points, as in the calculation of the normal state variational 
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Fig. 11. Energy gain per site AEg/Ns in the SDW state against electron density p is plotted for t' = (open 
diamond) and for t' ^ —0.10 (closed diamond). Error bar is about 0.0003. AEg/Ns in the superconducting 
state is also plotted again with the same symbols as in Fig. 9 with the diminished vertical scale. 



4-2 Phase Boundary between SC and SDW States 

As shown by the open diamonds in Fig. 11, the energy gain per site in the SDW state rises very 
sharply from p ~ 0.84 in the case of t' — 0. Already at p = 0.84 it takes 0.0023 and is slightly 
larger than that in the SC state. However at p = 0.80 there is no more stable SDW state. By 
extrapolating the sharply rising curve of the SDW energy gain per site as a function of p, the 
boundary between the SDW and the SC states is given at p = 0.838. Optimized values of M are 
plotted in Fig. 12. 

Since finite t' deteriorates the nesting property of the Fermi surface, the energy gain per site in the 
SDW state should decrease with introduction of negative t' . The calculated values for t' — —0.10 
is shown in Fig. 11 by closed diamonds. The energy gain per site in the SC state given by closed 
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Fig. 12. Optimized value of SDW gap parameter M is plotted as a function of p for two values of t' = 
(open diamond) and —0.10 (closed diamond), corresponding to two SDW curves in Fig. 11. 



circles increases when t' = —0.10. The phase boundary between the two is estimated at about 
0.873 as the intersection of the extrapolation of the line linking the two closed diamonds and the 
curve for the SC energy gain per site with t' = —0.10. So the boundary is pushed up by Ap = 0.035 
but not enough as to move the boundary up to p = 0.95, the observed boundary in La2-xSr3;.Cu04. 
The obtained phase diagram in the t'-p plane is shown in Fig. 13. 

§5. Comparison with Experiments and Other Results 

5. 1 Energy Gain per Site in the SC State 
By means of the cell-perturbation method, Feiner et al. were able to reduced the d-p model into 



an effective one-band Hubbard model |2^. In their result the transfer energy t slightly depends 
on the occupancy on the two related sites. Using typical parameter values for the d-p model t is 
estimated at 0.51 eV for hole transfer. U is estimated at 3.4 eV so that U/t = 6.7. With this 
relatively small U the authors derive an appropriate value of exchange interaction constant J for 
the t-J model, appearing later, thanks to the ferromagnetic interaction term and the Coulomb 
interaction between nearest neighbor sites in the effective one-band Hubbard model. These values 
of t and U are common for the cuprate high-T^ superconductors. The t' value was estimated at 
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Fig. 13. Phase diagram of superconducting (SC) and SDW phases in cuprates high- Tc superconductors in 
the plane of temperature T versus electron and hole densities. The thick curve was decided as written in 
the text. The dotted line is an extrapolation. 



around —0.06 eV, i.e., —0.12 in units of t, for hole-type superconductors. They observed that its 
weak dependence on each cuprate is correlated with the values of Tc among the cuprates. Since 
the additional spin-dependent ferromagnetic term and nearest neighbor Coulomb interactions are 
relatively small, the 2D Hubbard model without these additional interaction but with the above- 
mentioned parameter values should give a right magnitude of the SC energy gain if it actually gives 
a finite value. 

For [/ = 8 and t' = the maximum SC energy gain per site is approximately given by the value 
at p = 0.84 in the neighborhood of the SC-SDW phase boundary. It is about 0.0010 as the average 
of the values for both shell states. For U ^ 6.7 this value should diminish by a factor 3/4 according 
to Fig. 4 so that the calculated maximum SC energy gain per site is 0.00038 eV/site, if we take 
t = 0.51 eV. The value of t' for hole conduction was evaluated by several authors at, e.g., —0.06 eV 
2^,^, —0.124 eV (for the La system) |]26|, and —0.17 eV for hole-doped cuprates. If we take 



t' = —0.06 eV = —0.12 given by Feiner et al. the above energy gain should be multiplied by a 
factor 2.0 and becomes 0.00076 eV/site. Incidentally, this value of t' is in the range of the optimal 
energy gain, t' = —0.1 ~ —0.15. 
According to Hao at al. the critical magnetic field Hc{0) at zero temperature was estimated to 
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be 1.10 Tesla for YBa2Cu307 [2S]. From the expression H'^/'^ti for the condensation energy in the 
SC state, the experimental value for the SC energy gain per site is equal to 0.00026 eV per Cu 
site in layers. Triscone et al. gave slightly larger values of for two samples of YBa2Cu307, 
1.231 and 1.364 Tesla, corresponding to condensation energies 0.00033 and 0.00040 eV/(Cu site), 
respectively. 

Another source of information on the condensation energy is the specific heat reported by Loram 
et al. on YBa2Cu307 By numerically integrating the SC specific heat minus the normal-state 
one with respect to temperature from zero to just T^, the condensation energy was obtained as 
0.000168 eV/(Cu site). Since the appreciable fluctuation contribution is observed at temperatures 
above T^, which is recognized as the maximum position of the specific heat, and since it was simply 
neglected, this value is in fair agreement with the value given from the critical field. 

The above-mentioned calculated value of the maximum SC energy gain per site is in reasonable 
agreement with the experimental SC condensation energy in view of simplifications and uncertain- 
ties in the parameter values of the model. This agreement strongly indicates that the 2D Hubbard 
model includes essential ingredients for the SC in the cuprate superconductors. 

The present results of the energy gain per site in the SC state can be compared with the results 
for the t-J model. For the value of exchange interaction constant J = 4t^/f/ = 0.5 corresponding 
to [/ = 8, Fig. 7 in |^ by Yokoyama and Ogata allows to estimate the SC energy gain per site 
calculated for the t-J model. At p = 0.84 it is 0.026t. This is 17 times larger than that obtained for 
the 2D Hubbard model at the same electron density in §3. If we take 0.0010 for the SC energy gain 
per site of the 2D Hubbard model, judging from the average of the open and closed-shell curves, 
the ratio is 26. Since the 2D Hubbard model has been found to give a sound SC condensation 
energy as seen above, we have to judge that the fault is on the 2D t-J model as an effective 
Hamiltonian of the 2D Hubbard model, i.e., it gives too large an SC condensation energy at least 
in the parameter region where J/t ~ 1/2. This means that the t-J model made of the leading two 
terms in the expansion in powers of t/U of the canonical transformation of the Hubbard model 
should be treated together with the higher-order terms for it to give a realistic SC condensation 
energy. 

The t-J model is derived also starting from the d-p model or three-band Hubbard model. As 
such its parameter values are derived in literatures. The values are given in wide ranges, even by 
a single group. The value of t is estimated at 0.224 eV by Tohyama et al. for the La system [p6| . 



at 0.51 eV by Feiner et al. in the range of 0.29 ~ 0.98 eV by Batista et al. [|2[- The value of 



J is obtained at 0.13 eV for La2Cu04 from Raman scattering data |3^. Theoretical values given 
by these authors for typical cases are around this value. If we take J = 0.13 eV and t = 2 J, the 
SC energy gain of the 2D t-J model with p = 0.84 is estimated at 0.026t = 0.0068eV/site by using 



31I]; it is 26 times larger than the value obtained from Hao et al.'s He value. Around this value of 
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J /t = 0.5 the SC energy gain of the t-J model is considered to be roughly proportional tot-J/t = J, 
if the SC energy gain of the model with normalized coupling constant J/t is roughly proportional 
to J/t. Then, the above-mentioned large ratio applies to the most parameter sets. Only when J/t 
becomes so small, e.g., with increase of t, that SC almost vanishes, the above ratio may decrease 
to unity. Thus with plausible parameter sets, the t-J model is very probable to give too large a 
SC condensation energy. Again the higher order correction terms must not be neglected, for the 
model to work properly at least concerning the SC condensation energy, so that they oppose to the 
occurrence to SC P^; the present result does not support the correction terms of the type that 



facilitate SC 34 



5.2 Phase Diagram in the Tc-versus-p Plane 

According to our results in §4, in the case of small negative t' the SC region extends from p ~ 0.76 
to p ~ 0.86 with smoothly increasing SC condensation energy with increase of p. This suggests 
that in this p region Tc is finite and smoothly increases with increase of p. The range of this p 
region and the expected smooth increase of are in fair agreement with the features of the T^- 
versus-p {p = 1 — p, doped-hole density) phase diagram in the optimal to overdoped region, i.e., 
0.15 < p < 0.25 or 0.75 < p < 0.85, where Tc increases with increasing p [p3 |. 



Corresponding to the underdoped region with 0.05 < p < 0.15, the present results do not give an 
SC region but give an SDW region. However, there is a possibility that more elaborate calculations 
give a kind of SC phase where SC and SDW coexist and that the phase diagram in the underdoped 
region is also provided by the 2D Hubbard model. The reasons are as follows: Giamarchi et al. 
showed that a uniform wave function describing the coexistence of SC and SDW gives a lower 
energy than the pure SDW wave function in the system of 32 electrons on the 6x6 lattice 0. We 
have checked on the 10 x 10 lattice that such a coexistence wave function is possible only in the p 
region higher than the above-mentioned boundary between SC and SDW. Therefore, in the higher 
p region there may be a coexistence phase in a certain region. However, there are no experiments 
confirming a uniform coexistence except for recent reports suggesting a non-uniform coexistence 
in Lai.6-xNdo.4Sr^Cu04 [Q. Further, stripe-type solutions are known to give a lower total energy 
for the SDW state in the low hole-doping region Therefore, the above-mentioned uniform 

coexistence state may be dominated by a more complicated non-uniform way of coexistence of 
SC and SDW. Such a ground state could be obtained only if we use sufficiently large sizes of the 
system, as is the case with the SDW stripe solutions ||36|. With this speculated coexistent SC 



state, we expect that is maximum at the above-mentioned phase boundary at p ~ 0.86, since 
the plausible coexistence state should have lower than the value of at the boundary due to a 
smaller electron density carrying SC properties. At temperatures higher than Tc, the SDW ordering 
without the SC one may remain. Such a coexistent SC state may be consistent with the feature of 



the Tc-p phase diagram in the underdoped region [23|. Of course, all these possibilities remain to 
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be confirmed with actual calculations. 



5.3 Comments on the Quantum Monte Carlo Results 

Whether the 2D Hubbard model drives SC or not has been an important unsettled issue for 
theoretical and computational physics. Although quantum M. C. studies recognized the enhance- 
ment of the SC susceptibility, dominance of SC has not been confirmed. However, the numerical 
studies have been under severe restrictions to the parameter values in which reliable results were 
obtainable, e.g., with respect to the system size, value of [//t, temperature and so on. Within these 
few years the occurrence of the SC was strongly indicated when t' is introduced |^,^|. When the 
shell structure approaches to the open one, even if it was closed, SC correlations were observed to 
be much enhanced On the contrary, recently Zhang et al. showed that SC features weaken 
with increase of the system size with U up to 8t |3^. However, the state with p = 0.85 and t' = 
with which their main results are exhibited is located in the SDW region according to our results. 
With this value of p the dominant SDW is considered to have diminished the SC correlations for 
some values of U. Even in the case where the coexistence of SC and SDW might have occurred, the 
SC features should have been different from what we expect in the usual uniform superconducting 
state. Therefore, their results do not disprove the existence of SC in the p region where we find SC 
dominating. Further the quantum M.C. calculations were done for closed-shell electron densities. 
This shell state must have led to a very small energy gain and weakened the SC features when the 
system size is not sufficiently large. The distribution of the single particle energy ej^ around the 
Fermi energy are scarce in the cases they treated. These reasons are suspected to have made the 
SC features in their work very weak. A study with an elaborate method allowing high precision is 



under way searching enhanced SC correlations with appropriate values of p 
§6. Summary 

Using the variational Monte Carlo method we have investigated the dependencies of the energy 
gain per site in the SC state on important parameters such as on-site Coulomb energy U, next 
nearest neighbor transfer energy t' and electron density per site p. We worked mainly with the 
10 X 10 lattice. The energy gain is maximized at about U = 8. It is maximized for t' = —0.10 ~ 
—0.15, reaching twice the value at t' = 0. This increase of the SC energy gain was ascribed to the 
fact that a high density of the one-electron levels around the van Hove singularity go close to ep due 
to an appropriate negative value of t'. The SC pair correlation functions increased correspondingly. 
With U = 8 and t' = —0.10 or —0.15, the energy gain starts to be finite at p ~ 0.68, smoothly 
increasing with increase of p from p = 0.76 up to p = 0.92 up to which we have made calculation. 
The curves of the energy gain as a function of p have only a weak dependence on the shell structure 
or on whether some of the fc-points are partially filled or not (open or closed shell) in the U = 
limit electron distribution. This definitely suggests that the curves are already close to the bulk 
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limit ones. On the other hand, with U — 8 and t' — 0, the energy gain in the open shell state 
is sizable but that in the closed shell state is rather suppressed so that the obtained values seem 
to have still appreciable system size dependence, although the average curve between the open- 
and closed-shell curves is considered to be a fair approximation to the bulk hmit. We have also 
calculated the energy gain in the commensurate SDW state. It sharply rises as a function of p 
starting at about p ~ 0.84 for p — 0, making the SDW state the lower energy state in the higher p 
region above this value. A finite value oit' ~ —0.1 brought about an upward shift of the boundary 
by At' ~ 0.035. Therefore, our calculations confirmed that the electronic repulsive interaction can 
drive SC by itself in the region from p ~ 0.76 up to the boundary at p ~ 0.87 with t' ~ —0.1. The 
maximum SC energy gain per site in the SC region was found to be close to the experimental SC 
condensation energy evaluated from He and the specific heat of YBa2Cu307. On the other hand, 
the corresponding t-J model was pointed out to give a value larger enormously by a factor exceeding 
20, indicating that it is not quantitatively reliable as a model for high-Tc cuprates without taking 
account of higher-order correction terms. The feature of the experimental phase diagram in the 
Tc-p plane in the region of p lower than the p of the maximum Tc was asserted to be in a good 
agreement with our results. These results strongly suggest that the simple 2D Hubbard model 
includes essential ingredients of high-Tc cuprates. A speculative argument on the higher p region 
outside our SC region was given concerning the coexistence of SC and SDW. 
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